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skRODUCnON 

"  Ultxasomd  B-scan  imaging  is  now  a  well  established  and  valuable  clinical 
tool.  Inprovements  in  transducer  arrays  and  microprocessor  controls  have  lead 
to  the  development  of  real  time  linear  and  sector  scanners  which  produce  images 
of  remarkable  clarity  and  resolution  conpared  with  B-scanners  of  only  a  few 
years  ago.  Further  inprovements  in  B-scan  images  are  predicted  to  occur  as 
larger  transducers  apertures  and  improved  dynamic  focusing  methods  are  enployed. 

The  use  of  Doppler  xiltrasound  alone  or  in  combination  with  real  time  B-scan 
imaging  is  expected  to  increase  in  inportance  as  the  clinical  significance  of 
hi^  resolution  Doppler  images  is  appreciated. 

During  the  past  decade  X-ray  C.T.  (conputed  tomography)  and  recently  NMR  imaging 
have  made  remarkable  contributions  to  the  field  of  diagnostic  imaging.  X-ray  (C.T.) 
and  NMR  images  not  only  provide  resolution  of  about  1  nm  but 

also  can  be  calibrated  by  absolute  reference  standards.  The  resulting  quantitative 
images  have  proven  to  have  valxjable  diagnostic  valxie  because  of  the  greater 
ability  thus  provided  to  distinguish  healthy  and  diseased  tissue  by  their  image 
values .  ^ 

It  is  nat\iral  to  ask  if  the  successful  B-scanner  technology  could  be 
improved  further  by  incorporating  quantitative  tissue  characterization  features . 

Such  tissvie  characterization  features,  based  on  back  scatter,  are  now  being 
investigated  for  inclusion  on  clinical  B-scanners,  but  the  tissxje  characterization 
they  provide  is  based  on  structural  and  statistical  properties  of  tissxaes  and  not 
absolute  tissue  properties.  The  statistical  properties  of  tissues,  e.g.,  the 
spatial  Fourier  transform,  are  often  correlated  with  the  state  of  health  or  disease 
and  are  therefore  valuable  but  they  are  not  easy  to  measure  quantitatively  with 
back  scatter  techniques  mless  special  precautions  are  taken  to  model  or  compensate 
for  angular  dependence.  These  back  scatter  statistical  properties  are  dependent 
on  the  gradient  of  the  density  of  the  tissijes  and  thus  the  constant  of  integration 
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has  been  lost  [  1  ] .  Thvis  we  have  partially  answered  our  question.  We  have  noted  that 
sone  quantitative  information  about  tissvie  properties  can  be  obtained  from 
B-scans  but  not  the  absolute  mechanical  properties  of  tissue.  It  is  then  appropriate 
to  ask  if  ultrasomd  can  in  principle  provide  additional  tissue  properties.  The 
answer  is  yes  if  the  scattering  into  many  angles  is  used.  We  next  describe  how 
such  scattering  has  been  investigated  for  imaging  applications. 

During  the  past  decade  many  investigators  have  attenpted  to  develop  an  tiltra- 
somd  tomography  vdiich  could  provide  resolution  and  absolute  tissue  characterization 
analogpus  to  x-ray  C.T.  The  first  ultrasound  C.T.  systems  appeared  at 

aboxit  the  same  time  that  NMR  imaging  was  demonstrated  and  shortly  after  x-ray  C.T. 
was  announced  [2,3].  The  development  of  tiltrasomd  C.T.  has  not  shown  the  rapid  develop 
ment  into  a  clinical  tool  as  have  the  other  modes.  The  rapid  development  of  IIMR 
imaging  has  been  siarprising  even  to  many  closely  associated  with  medical  imaging. 

The  reasons  vdiy  ultrasound  imaging  has  not  followed  the  same  pace  is  now  becoming 
more  clear. 

Ultrasound  C.T.  has  not  made  the  rapid  progress  of  x-ray  C.T.  or  NMR  imaging 
for  several  reasons:  (1)  the  theory  for  x-ray  C.T.  or  NMR  imaging  is  much  siapler 
than  that  of  ultreisomd  C.T. .  (2)  neither  x-ray  or  NMR  imaging  are  subject  to 
diffraction  and  refraction  effects,  but  these  two  effects  dominate  the  ultrasomd 
sigial  detection  process.  (3)  x-ray  detector  arrays  and  NMR  detectors  are 
relatively  sinple  to  design  and  build  compared  to  ultrasomd  detector  arrays. 

This  chapter  will  address  sioggestions  for  inproving  the  resolution  and  quantitative 
imaging  of  ultrasomd  C.T.  by  addressing  the  theory  for  iiltrasomd  C.T.  and  treating 
refraction  and  diffraction  effects. 

Initial  attenpts  to  develop  ultrasomd  C.T.  greatly  siaplified  the  problem 
by  treating  the  propagation  of  ultrasomd  as  if  the  acoustic  energy  followed 
straigjht  line  rays  like  x-rays.  As  a  conseqxience ,  refraction  and  diffraction 
effects  limited  resolution  to  about  10  to  20  wavelengths.  Correcting  for  refraction 
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of  the  ray  paths  did  not  help  because  resolution  is  determined  nostly  by  diffraction 
effects .  As  the  diffraction  problem  became  better  understood  several  investigators 
sxjggested  including  diffraction  and  refraction  in  the  ultrasound  imaging  formulations 
by  borrowing  Wolf's  clever  perturbatiai  theory,  inverse  Bom,  scattering  algorithm 
from  optical  holography  theory  [4,5] .  Mjch  effort  has  been  expended  in 
developing  this  approach  with  the  greatest  siaccess  coming  from  replacing  the  Bom 
approximation  with  the  more  accurate  ifytov  approximations.  Although  this  work  has 
brought  fresh  and  iiqportant  insight  to  the  quest,  there  are  several  serious  defects 
in  these  Rytov  or  Bom  perturbation  methods  [  6  ] :  (1)  these  methods  must  assime 
unrealistically  low  absorption.  (2)  the  refraction  is  only  treated  approximately, 
e.g.,  the  far  field  of  a  siaple  lens  is  not  given  correctly  by  either  the  Bom  or 
Rytov  approximate  methods.  (3)  the  inclusion  of  density  scattering  has  been  solved 
but  still  depends  on  attenioation  being  treated  correctly. 

The  use  of  an  attenuating  bath  has  been  suggested  to  allow  imaging  of  himan 
tissue  at  hi^er  frequencies,  e.g.,  1  to  3  mhz.  A  new  algorithm  to  allow  better 
treatment  of  refraction  effects  for  ultrasound  C.T.  has  been  proposed  [7  ].  Such 
approaches  may  help,  yet  the  fact  remains  that  the  perturbation  methods  are  not 
exact  except  in  the  limit  of  vanishingly  small  perturbations.  Therefore,  an  ultra¬ 
sound  imaging  method  based  upon  a  more  exact  treatment  of  the  wave  equation  would 
have  the  capability  of  imaging  objects  with  larger  attenua¬ 
tion  and  larger  variations  in  refractive  index  and  density.  Such  a  method  would 
also  have  the  desirable  feature  of  acting  as  a  standard  of  conparison. 

We  are  aware  of  tvro  methods  which  do  not  depend  on  the  object  producing  only 
small  perturbations  to  the  total  field.  The  first  method  is  based  upon  extra¬ 
polating  the  scattered  data  to  zero  wavelength.  This  method  looks  promising  for 
medical  imaging  and  has  potential  for  wider  development  and  use  if  its  sensitivity 
to  noise  can  be  reduced  [8,9]. 


The  second  method  is  less  sensitive  to  noise  and  provides  "exact"  inverse 
scattering  solutions  (by  exact  we  mean  the  equation  or  theoretical  nodel  used 
for  imaging  contains  no  unrealistic  or  major  approximations) .  The  second  method 
although  "exact"  in  the  above  sense  is  nevertheless  iterative  in  nature  and  may 
probably  reqiiire  a  later  generation  of  cocnputers  for  cost  effective  inplementation. 
Nevertheless,  the  method  is  of  academic  and  laboratory  significance  and  may  serve 
as  the  starting  point  for  further  inprovements  in  \altrasonic  imaging;  a  description  follows 


ULTRASONIC  imCING  BY  SOLUTION  OF  THE  INVERSE  SCATTERING  PROBLEM 


The  ultrasound  imaging  problem  in  its  most  general  and  complete  form  is 
equivalent  to  solving  the  so-called  inverse  scattering  problem.  The  inverse  scat¬ 
tering  problem  may  be  defined  as  finding  the  spatial  distribution  y  of  material 
properties  so  that  the  wave  eqijation  scattered  field  solutions  f ^  ^  for  any  incident 
field  f^^^  matches  the  given  scattered  field.  The  solution  to  the  direct  scattering  problem 
is  defined  to  be  the  solution  f '  '  (x)  to  a  wave  equation  at  some  point  x  for 
y  given.  Thus,  the  solution  y  to  the  inverse  scattering  problem  is  that  y 
for  \i^ich  the  solution  to  the  direct  scattering  problem  matches  given  scattered 
field  valvies.  In  the  usual  ultrasound  imaging  environment  only  the  scattered 
fields  on  a  detector  are  given  by  measurement,  but  the  scattered  fields  within  the 
object  are  not  known.  If  the  fields  were  knovn  within  the  object  then  the  solution 
to  the  inverse  scattering  problem  could  be  obtained  by  a  linear  operation.  Since 
both  the  internal  field  and  material  property  are  unknovm,  the  solution  to  the 
inverse  scattering  problem  also  reqiiires,  at  least  inplicitly,  the  solution  to  the 
direct  scattering  problem. 

Because  of  the  past  difficvilty  of  solving  the  inverse  scattering  problem, 
approximate  solutions  have  been  developed.  (>ie  of  the  most  widely  known  approximation 
methods  is  based  upon  using  first  order  perturbation  theory  to  linearize  the  procedure 
for  solving  the  wave  eqxjation  as  a  function  of  y  .  This  approach  leads  to  the 
derivaticxi  of  the  well  known  Bom  and  Rytov  approximations  [5].  For  the  Bom 


approximation,  eqtiations  for  the  perturbation  in  the  acovistic  pressure  field  p 
and  perturbation  in  material  properties  y  are  given  by  the  equations 


(1) 

P  *  Po  Pi 

CM 

Y  -  Yq  + 

viiere  p  is  the  total  field  when  Y  ®  Yq  +  and  p2  is  the  field  vdien  Y  =  Yq 
vdiere  Yq  is  a  constant  and  is  the  material  property  of  the  homogeneous  background 
I  5  ] .  The  Rytov  approximation  is  obtained  by  the  transformations 

(3)  p  =  exp(W) 

(4)  W  =  Wq+Wj^ 

viiere  Wq  is  the  logarithm  of  the  incident  field  when  y-^  =  0  ,  where  W  is  the 
logarithn  of  the  total  field  vdien  Y  “  Yq  +  Y]^  and  where  is  the  pertxjbation  in 
W  caused  by  the  perturbation  Yj^  in  y  • 

For  medical  applications,  the  assui5>tion  that  Yj^  ,  Pj^  ,  and  W  are  small  with 
respect  to  Yq  ,  Pq  f  and  Wq  ,  respectively,  may  in  many  cases  be  erroneous  because 
of  refraction  and  attenuation.  This  may  be  seen  by  inspection  of  the  Helmholtz  wave 
equation 

9  2 

(5)  p(x)  +  — j  p(x)  =  0 

“  lc(x)]^  - 

The  speed  of  sound  c(x)  may  be  considered  to  be  coqplex  and  we  could  write 

_2 

c(x)  «  iCj(x)  •  However  it  is  more  convenient  to  first  expand  [c(x)]  by 
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equation  (2)  before  examining  the  real  and  imaginary  parts.  By  adding 
2  2-2  -2 

kQF(x)p(x)=u)  (Cq  -[c(x)1  )p(x)  to  both  sides  of  equation  (5)  we  obtain 

(6)  p(x)  +  k2[l  -  F(x) 3p(x)  =  0 

2  2  2 

v^re  kQ  *  0)  /Cq  is  the  background  wave  nunber  sqxjared.  Thus  by  inspection  of 

2  2 

equation  (2)  we  see  that  Yq  “  kg  and  =  F(5.)kQ  .  We  now  consider  two  exanples 
vMch  illustrate  the  effect  of  F(x)  on  the  relative  sizes  of  pg  and  p^  . 

The  first  exanple  is  that  of  a  siirq>le  lens.  It  is  clear  that  a  small  value 
of  F(x)  ,  say  five  or  ten  percent,  distributed  in  a  lens  shaped  region  can  focus 
an  incident  plane  wave  to  a  local  zone  where  the  total  field  p  may  be  many  times 
larger  than  the  ur^erturbed  field  pg  .  Thus,  the  scattered  or  perturbed  field 
Pj^  can  not  be  small  in  conparison  to  Pg  .  The  second  example  is  that  of  a  plane 
vave  normally  incident  upon  a  block  of  homogeneous  material  with  a  small  but  finite 
loss  corresponding  to  the  imaginary  part  of  F(x)  .  It  is  easily  shawn  that  even 
thoui^  the  imaginary  part  of  F(x)  is  small  the  line  integral  of  the  linear 
acoustical  attenuation  coefficient  can  reduce  the  anplitude  of  the  internal  field 
to  a  small  fraction  of  the  incident  field  in  only  10  cm.  of  travel.  Soft  tissue 
has  an  attenuation  of  only  about  1  db  MHz  ^  .  For  an  average  speed  of  sound 

of  1500  m.  sec  ^  and  at  1  MHz,  1  dbMHi?'  cm”*^  corresponds  to  a  ratio  of  imaginary 
to  real  part  of  the  wave  number  of  about  0.005.  Nevertheless,  at  1  MHz  at  a  depth 
of  10  cm  in  the  block  the  field  is  down  10  db  relative  to  the  field  at  the  surface, 
that  is,  the  pressure  is  reduced  by  68  percent. 

These  exanples  illustrate  the  need  for  caution  when  applying  the  Bom  or 
Bytov  theory.  The  situation  can  be  improved  if  the  background  fluid  is  tailored 
to  have  a  real  and  imaginary  speed  of  sound  equal  to  the  mean  respectively  of  the 
real  and  imaginary  parts  of  the  speed  of  sound  for  the  object  in  the  background 
fluiid.  A  method  for  obtaining  an  inverse  scattering  solution  to  the  exact  Helmholtz 


wave  equation  or  vAiich  does  not  require  act\jal  field  to  be  nearly  equal  to  the 
incident  field  with  no  object  present,  wauld  avoid  difficulties  illustrated  by  the 
two  exanples  above.  We  next  review  the  background  of  a  method  for  obtaining  a 
direct  solution  to  the  exact  Helmholtz  wave  equation  \diich  we  will  use  to  develop 
a  robust  inverse  scattering  solution  to  the  exact  Helitiholtz  wave  equation. 

Methods  for  solving  the  direct  scattering  problem  by  the  method  of  moments, 
of  vdiich  the  Galerkin  method  is  a  special  case,  are  presented  in  the  book  by 
Harrington  [10]  and  the  paper  by  Richmond  [11].  Recently  the  method  of  Richmond 
was  adapted  independently  by  Hagnann  [12]  and  by  Yoon  et  al  [13]  to  solve  the 
inverse  scattering  problem  for  the  case  of  a  single  source.  They  found  that  the 
inverse  scattering  problem  for  a  single  source  is  an  ill-posed  problem  and  therefore 
very  small  amomts  of  noise  in  the  data  caused  a  large  amomt  of  noise  in  the 
solution  for  y  .  As  the  nunber  of  points  in  the  image  domain  grows,  the  condition 
nunber  of  the  matrices  to  be  inverted  by  the  single  source  method  grows  rapidly  and 
the  solution  becomes  unstable  and  inacctarate.  The  stability  of  the  solution  also 
decreased  as  the  detector  to  object  distance  was  increased. 

The  ill-posed  nature  of  the  inverse  scattering  solution  may  be  largely  rsnoved 
by  solving  a  system  of  nonlinear  equations  resulting  from  using  multiple  sources  [14] . 
The  multiple  source  solution  has  demonstrated  good  tolerance  to  additive  noise  in 
the  scattering  date  and  is  not  sensitive  to  the  position  of  the  detection 
for  over- determined  cases. 

FORMULATiai  OF  EQUATIONS  WHICH  MAY  BE  SOLVED  FOR  DIRECT  AND  INVERSE  SCATTERING  SOLUnONS 

We  now  introduce  an  inproved  wave  equation  containing  equation  (6)  as  a 
special  case  and  introduce  an  equivalent  integral  equation.  We  then  transform  the 
integral  equation  into  a  set  of  multivariate  quadratic  equations.  The  solution 
to  these  quadratic  equations  lends  to  the  direct  scattering  solution  and  the 


inverse  scattering  solution. 
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The  wave  equation  for  a  body  iimersed  in  a  homogeneous  fluid  is  given  by 


I  1] 

(7)  +  k^f^Cx)  =  (x)f^(x)  -  V*[Yp(x)Vfp(x)] 


where  f^Cx)  is  a  scalar  field  such  as  acoxistic  pressure  for  a  source  at  position 
0  ,  Icq  *=  u/cq  ,  ui  =  angular  frequency  of  incident  field,  Cq  is  constant  speed  of 
sound  in  a  homogeneoiJis  background  fluid,  =  (<~  ~ 

the  fractional  change  in  conppessibility  k  and  density  p  ,  respectively  outside 
the  body  in  the  fluid  background.  A  useful  integral  equation  vdiich  is  equivalent  to 
(7)  and  includes  the  boundary  conditions  is  [  1  ] 


(8)  f^®\x)  ^  f0(2i) 


f^^^(x)  = 


[k^Y^  (x ' )  f 0  (x ' )  g  (x ' ;  x)  +  YpV’f 0  (x ' )  •  Vg  (x ' ;  x)  ]  d^  X ' 


viiere  x=  (Xj^ . x^)  ,  x'  •=  (Xp...,XQ)  are  spatial  coordinates  in  a  space  of 

dimension  Q  ^  3  ,  f^®^  (x)  is  the  scattered  field  f^Cx)  -  f^^^  (x)  ,  f^^^  (x)  is 

the  field  fp(x)  for  no  object  present;  i.e.,  for  Y^  =  Yp  =  0  ,  and  where 

(2) 

g(x' ;x)  is  a  Green's  function.  For  two-dimensional  problems  g(x' ;x)  =  Hq  (kglx-x 
where  is  the  first  order  Hankel  function  of  the  second  kind.  For  three- 

dimensional  problems  g(x' :x)  “  (^■'’’Ix-x' |)”^  e:qp(ikQ |x-x' |)  ,  and  \diere  1*1  means 
absolute  valiie. 

We  now  expand  the  fields  and  f^®^  and  factors  containing  the 

material  properties  using  an  accurate  basis  set  {i(;^(x) }  .  Let 


N 

(9,10)  fp(x)  =  a^.  i(;.(x) 


*>01  ♦i® 


(11.12) 


(13)  Yp(x')Vf0(x')  »  Aqp.  il>.(x’)iq  , 

where  N  is  the  nuiiber  of  basis  functicxis  and  Q  is  the  dimension  of  the  space 
and  corresponding  Greens  function,  and  vdiere  is  a  mit  vector  along  coordinate 
q  .  Upon  substitution  of  (13)  into  (8)  we  obtain  upon  interchaning  integration  and 
sunnation, 

(V  =03  Ji  Ji  V 


N 


N 


^ere 


(15,16)  gj(x)  =  Jg(x';x)’^j(x')d^x' 


G  .  (x)  = 

q]  - 


x’ 


3x_ 


Using  the  "expansion  method"  equation  (14)  coijld  be  developed  into  a  system  of 
algebraic  equations  describing  the  scattering  by  evalxjating  (14)  at  a  set  of  grid 
points  {x. }  .  Alternately,  by  taking  the  inner  product  of  (14)  with  all  members 
of  a  weighting  or  testing  set  of  basis  functions  {W^(x) }  a  system  of  N  algebraic 
equations  may  be  generated  by  the  '*mBthod  of  moments"  [10].  If  V^(x)  =  i|^^(x)  for 
all  i  and  {ip^}  are  independent  then  the  procedure  becomes  the  classical  "Galerkin 
method"  [10].  If  W^(x)  =  6(x-Xj^)  the  method  of  moments  becomes  the  expansion 
method.  The  choice  of  the  expansion  basis  functions  {^^(x) }  and  the  testing 
functions  {Wj^(x) )  is  problem  dependent  and  often  significantly  effects  the  accxiracy 
of  solution  and  the  nixnber  of  e:qpansion  basis  functions  needed  to  obtain  a  solution 
[  10] .  The  sine  basis  functions  have  many  properties  vdiich  make  them  especially 
attractive  for  applications  using  the  method  of  moments  and  in  the  remaining  sections 
we  will  use  only  basis  sets  based  on  the  sine  functions. 
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ALGEBRAIC  SCATTERING  EQUATIONS  raj^VED  USING  SINC  BASIS  FUNCTIONS 

The  sine  functions  are  attractive  as  a  basis  set  because  v^en  properly 
formulated  the  basis  set  is  conplete,  orthogonal  and  retains  sane  local  or  regional 
properties  (by  local  or  regional  we  mean  that  each  coefficient  of  the  basis 
expansion  are  related  to  the  values  of  the  function  to  be  ejqjanded  in  a  unique  local 
region  and  this  property  often  aids  in  incorporating  measvtrements  or  boundary 
conditions)  for  representing  spatially  band  limited  functions. 

In  one-dimensions  a  function  f(x)  may  be  expanded  in  a  basis  set  consisting 
of  sine  functions  shifted  in  integral  multiples  of  step  size  h  as  follows 

00 

(17)  f(x)  C(f.h)  ^  I  f(«.h)S(«,h)(x) 

£=_oo 

vhere  the  shifted  sine  functions  are  given  by 

(18)  S(l,h)(x)  4  . 

If  f(x)  is  band  limited  to  contain  no  Fourier  conponents  greater  than  -rr/h 
^len  f(x)  =  C(f,h)  .  The  orthogonality  of  {S(£,h)(x)}  is  expressed  by 


(19) 


C(f,h)(t)S(Jl.h)(t)dt 


h  f(£h) 


The  nth  derivative  of  function  C(f  ,h)  has  a  sine  e:q>ansion  given  by 


(20) 


T1  / 

c(f,h)(x)  =  hT^  I  (I  ajf  f(iih)}s(j.h)(x) 


j*=-00  I 


(21) 


5(n) 

hi 


-^S(£,l)  (X) 
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In  particialar 


(22) 


0  .  =  j 

^  j 


The  cardinal  function  C(f  ,h)  may  be  integrated  by  use  of  the  following 
formula 


(23) 


#00  00 

C(f.h)(t)dt  =  h  ^  f(£h) 

-00  £=-“ 


These  and  other  useful  properties  of  the  sine  function  for  solving  integral 
equations  and  for  other  applications  are  given  by  Stenger  [15  ] . 

We  take  products  of  Q  sine  functions  to  form  basis  sets  for  problems  in  Q 
dimensions.  We  illustrate  this  for  two  dimensions  by  setting 


(24) 


’^i  (m.n)  ^  ^^1^  S  (n .  h)  (X2) 


to  obtain 


(25)  f0(Xj^,X2)  =  II  f0(mh,nh)S(m,h)(x,)S(n,h)(x2) 

^  m  n 


Here  i  =  i(m,n)  means  i  is  an  index  v^ich  indexes  point  (m,n)  in  a  finite 
two  dimensional  array  of  grid  points.  Thus  by  (9-12),  (24,25)  and  by  property 
(17)  we  obtain  relationships  for  the  coefficients  a^^  ,  ,  c^^  as  follows 


(26.27)  =  £j,(n*..nh) 


^0i(m,n) 


f^^^  (mh,nh) 


(28,29)  d.,,  .  =  f^®\mh,nh)  ,  c-.,  .  =  Y,(inh,nh)f«(mh,nh) 


The  relationship  between  ,  Yp  .  is  obtained  by  noting  that 

YpVt'Vg  ■  Yp(x')^  Of/3Xq)  Og/Sx^)  .  We  examine  the  factor  Yp(x') (9f/3Xq) 
and  note  that  Of/3x^)  has  a  sine  esqjansion  via  application  of  (20,21)  to  (23). 
The  resiilt  is  the  product  of  Yp  and  a  sine  expansion  of  (df/dx^  .  We  next 
perform  a  sine  expansion  on  the  product  by  use  of  (17) .  These  nested  expansions 
are  simplified  by  noting  that  S(Jl,h)(x  »  rh)  »  (vdiere  6^^  =  0  if  and 

6j^  “1  if  i  =  r)  .  ^plying  these  steps  for  the  terms  Yp(3f/9Xj^)  and 
Yp(3f/9x2)  we  obtain  on  setting  x^  =  x  ,  X2  “  y 

(30a)  Y.(x.y)  -  h‘^  H  Y.(iTih.nh)  [[6^^  f0(ilh,nh)]S(m,h)  (x)S(n.h) (y) 

(30b)  Yp(x,y)  — -  =  h"^  J  J  Yp(nih,nh)  I[6^^  f«(mh,£h)]S(m.h)(x)S(n,h)(y)  . 

We  now  have  obtained  ejpressions  for  a^^  ,  b^^  ,  c^^  ,  d^^  and  in  terms 

of  fp(mh,nh)  ,  ,  Y^(nih,nh)  and  Yp(inh.nh)  as  per  (26-29)  and  (30).  These 

eiqsressions  coxild  be  substituted  into  (14)  to  obtain  a  set  of  equations  with 
coefficients  i{'^(x)  ,  gj(x)  and  ^^.(x)  to  be  evaluated.  Instead  of  evaluating 
the  coefficients  at  all  grid  points  (iih,nh)  we  use  the  Galerkin  method  to 
expand  gj(x)  and  G^(x)  terms  of  'l'^(x)  •  Let 

(31,32)  gj(x)  -  J  ij;^(x)  and  Gqj(x)  -  J  ij;^(x) 

viiere  i|'^(x)  is  given  by  (24).  Thus,  by  (25)  the  coefficients  g^^^  and 
are  given  by 

(33,34)  gj^  -  gj(2^)  and 
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On  sxjbstituticjn  of  (26-29),  (30)  and  (31,32)  into  (lA)  we  obtain  a  set  of 


equations  involving  only 
2 

f^(nh,nh)  ,  Uq  ,  Y^(inh,nh)  ,  Yp(nih,nh)  ,  ,  and  • 

ThxjB  on  factoring,  we  obtain 


(35) 


I  *1  «  Z  (f0i  -  -  n  +  Z  Aq0j 

J  9 


:(i) 


\diere  the  sm  over  represents  the  conplete  ejqjressions  for  (30) .  According 

to  the  Galerkin  method  we  may  eliminate  the  factor  by  taking  the  inner  product 
of  (35)  with  il/j  .  Let  us  introduce  the  natural  notation  f^(mh,nh)  =  n)  ^  ^ 

Y^(dh,rh)  S  „)  4  Y^(„ ,  YpCnii.nh)  4  Ypi(„,n)  4  , 

Sj(r,t)i(m,n)  "  Sji  “  g(r,t)(m,n)  ®qj (r,t)i(m,n)  “  ^qji  ^  q(r,t)(m,n)  ' 

vdiere  i  =  i(m,n)  . 

Then  by  (19)  and  (30)  we  obtain 


f(s) 

^0(in,n) 


(36) 


(i)  2 

^|Z)(m.n)  ■  ^0(m,n)  *  %  V(r,t)  ^0(r,t)  S(r,t)(m,n) 

^p(r,t)  I  Ha.t)  ^l(r,t)(m.n) 

^0(r,il)  ^2(r,t)(m,n)l^ 


+ 


Equation  (36)  is  the  conplete  equation  for  describing  the  direct  scattering 
problem  with  y^  and  Yp  known  and  the  inverse  scattering  problem  with  y^  and 
Yp  inknown.  Note  that  the  suns  over  r  and  t  are  convolutions  because  of  the 
definitions  of  the  Green' s  function  g(x'  ;x)  .  The  sun  over  2.  is  also  a  convolu¬ 
tion  as  seen  from  (20-22).  The  suns  over  r,t  and  I  may  be  rearranged  so  that 
all  coefficients  of  each  f^^^  are  sunned  and  grouped  as  a  total  coefficient. 
Such  a  form  is  useful  for  direct  scattering  problems  v^ere  the  y,.  and  y  are 

K  p 


i 


known.  This  form  is  given  by 


f(s)  4  f  _  s  y  Y  e  + 

^0(m,n)  ^0(m,n)  ^0(m.n)  ^*^0  ^K(r,t)  S(r,t)(m,n) 

I  ^p(Jl,t)i(£,t)(m,n)‘'"  \£  ^p(r,£)2(r,£)(m,n)'^^0(r,t) 


We  next  present  methods  for  obtaining  the  direct  and  inverse  scattering 
solutions  from  (36)  and  (37)  for  the  special  case  of  Yp(x')  ®  0  .  Even  with  this 
restriction  it  will  be  clear  that  the  general  case  of  Yp(x')  0  be  obtained 
by  extension  of  the  method .  This  conceptually  simple  extension  requires  additional 
programning  conplexity  and  additional  memory  for  storing  the  coefficients 
and  is  therefore  reserved  for  future  publications. 


MEIHODS  FOR  SOLVING  M)DEL  EQUATIONS  FPR  CASE  OF  Yp  =  0 

In  formulating  the  direct  and  inverse  scattering  solutions  for  the  case 

Yp  *  0  it  will  be  useful  to  use  the  single  index  notion  i  to  represent  a  point 

x^  with  i  =  i(m,n)  ,  i.e.  =  (x^,yp  =  (nih.nh)  since  the  more  explicate 

two  index  notation  (m,n)  is  now  not  required.  We  also  emphasize  the  major  role 

vdiich  subscript  0  will  now  play  in  the  inverse  scattering  solution  to  obtain 

unique  solutions.  In  order  to  obtain  a  sufficient  nunber  of  equations  to  solve 

for  f^®^  it  will  be  necessary  to  also  solve  for  f0(x)  inside  the  body.  Thus 

(s) 

we  must  consider  two  sets  of  equations,  namely,  one  set  relating  f^  anywhere 
to  f^  inside  the  body  and  a  second  set  relating  at  point  x^  in  the  body  to 
all  other  points  Xj  in  the  body.  It  will  also  be  useful  to  write  these  equations 
in  terms  of  residuals  and  r^^  .  Thus  we  write 

<3s> 


14 


1 


2nifi  ftn  '<3  0j  mj 


0-1 . * 


2  2 

\ihere  »  kg  and  »  kg  .  The  subscript  tn  in  (39)  refers  to  a 
measxireniBnt  position  on  a  detector.  The  subscript  i  in  (38)  refers  to  any  grid 
point  in  a  region  including  the  body.  By  inspection  a  necessary  condition  for 
existance  of  a  comnon  solution  to  (38)  and  (39)  in  that  the  nuiiber  of  equations 
exceed  the  nvniber  of  unknowns,  i.e. ,  that  M  ^  N  . 

In  his  stiidy  of  the  direct  scattering  problem,  Richnond  [11]  obtained  a 

similar  set  of  equations  to  (38)  and  (39)  above  using  pulse  basis  functions.  The 

advantage  of  the  sine  basic  functions  lies  in  their  economy.  Riohnond  recomnends 

10  sanple  points  per  incident  wavelength  for  his  evaluation  of  and  to 

be  accurate  for  problems  of  low  contrast,  i.e.,  |Y|^|  £  .10  .  We  find  that  4 

sanples  per  incident  wavelength  is  sufficient  for  accurate  recons  trtjct  ion  of 

for  objects  vdiose  spatial  frequency  spectnxn  of  does  not  exceed  the  reciprocal 

of  the  incident  wave  length.  For  solving  the  direct  scattering  problem,  i.e., 

finding  f^  ,  we  follow  the  steps  in  Richmond's  method,  namely:  first,  for 

given  solve  the  linear  system  (38)  for  f^^  ,  second,  for  f^^  from  step  one, 

(s) 

and  for  given  y^^  confute  f^*^  by  evaluating  the  sxm  in  (39) . 

Solving  the  inverse  scattering  problem  is  more  difficult  since  both  Yj  and 
fjjj  are  unknown  and  are  present  in  (38)  and  (39)  as  qmdratic  terms  or  products 
Yj  fg)j  .  We  have  investigated  [16 ,17  ]  three  methods  for  solving  this  nonlinear 
problem:  (1)  iteratively  solving  alternate  linear  systems  method,  i.e.,  solving 
of  (38)  for  fpj  with  y^j  fixed  then  solving  (39)  for  y^j  with  f^^  fixed, 

(2)  nonlinear  row  action  methods  [18  ] .  (3)  and  optimization  techniques  [19  ,20  ] . 

The  last  two  methods  are  conplex  and  have  not  been  cociputationally  faster  than  the 
first  method  [17  ]  so  we  only  describe  the  alternate  iterative  method  (or  AIM  method 
for  short). 


The  AIM  method  for  iteratively  solving  alternate  linear  systems  may  be 
described  by  the  following  algorithm:  Let  f^  or  y  be  the  kth  iterate  of 
f^  or  Y  respectively.  Then 

1.  Pick  a  trial  value  for  y^^^(x)  acciiracy  6  . 

2.  Solve  (38)  for  with  y^^^(x)  fixed. 

3.  Solve  (39)  for  y^^^^  with  f^^^\x)  fixed. 

4.  Test  for  convergence  by  forming  the  objective  function. 


(40) 


If  pQgj  <  6  ,  go  to  step  5,  else  go  to  step  2. 

5.  Stop.  Print  results,  make  file,  display  image,  etc. 

Various  methods  for  solving  the  linear  systems  may  be  employed  such  as  Gauss 
elimination,  fixed  point  methods  (applicable  only  to  equation  (38)),  e.g. 


(41) 


f(^l) 


N 


Oc) 


-.Uj 

3=1 


fOc) 

% 


or  the  Kacanarz  method  [21 ,22  ] . 


RESULTS  OF  CXJMPUTER  SIMULATION  STUDIES 

The  performance  of  the  5-step  algorithn  for  obtaining  y^  for  the  case 

Yp  =  0  was  investigated  using  an  11  by  11  pixel  test  object.  Simulated  scattering 

data  was  obtained  by  the  following  two  steps:  First  solve  equation  (38)  for  f^^^^  , 

with  Y^j  given  by  the  values  of  the  test  object.  Second  using  the  values  so 

(s) 

obtained  and  the  given  values,  evaluate  f^'  from  equation  (39).  The  values 

of  thus  obtained  represent  noise  free  data  or  f^^  .  The  distribution  is 

normalized  such  that  •  Af  '  l/"  1  “  ®  "‘"l  norm" 

(also  called  norm,  Ifllj^*  llf^D  •  Thus  noisy  data  f^®^  *  ^IR^  is 

obtained  for  algorithn  evaluation. 
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Figure  1  shows  the  11  by  11  test  object.  Figures  2-3  show  the  inverse  scattering 
soluti<xi  obtained  using  noisy  scattering  data  with  B  =  .10  ,  the  data  over- 
determined  such  that  *  2.AN  ,  and  the  inagp  band  limited  to  have  no  spatial 
frequencies  higher  than  the  reciprocal  of  the  incidait  fields  wave  length.  The 
detector  radius  was  10  incident  field  wave  lengths  and  the  pixel  separation  was 
one  fourth  incident  field  wave  lengths.  If  a  picture  quality  index  or  "one  norm 

error  rate"  a  =  *’'Ypicture“^TRUE*l^*^TRUE*l  ccnputed  it  is  seen  that  a  .06  . 
Without  the  spatial  filtering  a  .17.  Ihese  results  suggest  that  with  over- 
determination  of  the  data,  with  small  contrast  i.e.,  |y|  <_  .10  ,  with  a  small 

detector  radius ,  and  with  a  small  nuiiber  of  pixels  (here  N  =  11  x  H  =  121)  the 
problem  is  very  well  posed. 

The  question  of  the  effect  of  increasing  the  niniber  of  pixels  N  and  the 
contrast  in  y  on  the  noise  properties  of  the  solution  image  was  investigated. 

We  also  investigated  the  effect  of  increasing  the  detector  radius.  This  investiga¬ 
tion  was  conducted  to  determine  the  inverse  scattering  solution  when  the  internal 
fields  were  known  with  certainty.  An  image  computed  using  the  5-step  algorithm, 
could  not  have  less  noise  than  the  a  method  using  the  correct  internal  fields. 

The  solution  to  the  inverse  scattering  problem  with  the  correct  internal  fields 
given,  i.e. ,  evalviating  (39)  for  y  ,  we  call  the  "pseudo  inverse  scattering 
problem"  [23] .  The  sinxaltaneous  solution  of  equations  (38)  and  (39)  we  call  the 
"conplete  inverse  scattering  problem"  [23] . 

Our  stxjtfy  was  composed  of  the  following  6  steps  vdiich  provide  two  different 
upper  bounds  or  estimates  on  the  noise  present  in  the  solution  to  the  conplete 
inverse  scattering  problem. 

1.  Generate  the  scattered  field  data  vdiich  vrould  be  observed  at  a  detector 
using  a  (Saussian  object  distribution.  Define  f^®^  as  the  composite 
colinn  vector  of  <I>M  components  given  by 


(A2a) 


f(s)  4  f(s)  f(s)  f  (s)  i"] 

-  ^-1  ’  -2  »  •  •  •  •  i0  ^ 


\d[iere 


(42b) 


J<s)  4  (f(s=)(^) . . 


2.  Generate  a  Gavissian  randan  noise  vector  with  a  mean  of  zero 

and  normalized  such  that 


I. 


3.  Starting  with  the  overdetermined  linear  system  (39) 


Al  - 


the  least  squares  solution  for  x  '^ich  we  accept  as  "the  best 
solution,"  is  obtained  by  solving: 


A«Ax  -  AHf(*> 


Here  A"  is  the  conplex  conjugate  transpose  of  A  . 

4.  Compute  the  "upper  bound  on  the  norm  error  ratio  in  Ay  ,"  ,  using 

the  condition  nixtiber  inequality  [  23,  page  431] . 


II  Axil,  „  lA*^Af^®^ll,  . 

H  H 

Here  K2^(A  A)  is  the  one  norm  condition  number  of  A  A 
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(47) 


5. 


Determine  the  actual  norm  error  from  the  true  value  of  y  (x^ ) 
^miE  ’  solving  for  +  4^)  =  ^  from  the  equation 


i.e, , 


A^A(x,j^  +  Ax)  =  A”(f®  +  Af®) 


6.  Determine  Ax  by  solving  (47)  for  (Xtpuf.  ■*"  Ax)  compute  the  "one 
norm  error  ratio  of  the  image"  or  ct  by 


(48) 


‘^I"l  A 


From  this  type  of  analysis  of  the  pseudo-inverse  problem,  we  have  found 
the  following  trends. 

U 

1.  For  a  given  grid  spacing  h  ,  the  condition  nuiiber  Kj^(A  A)  is  an 

exponentially  increasing  function  of  the  edge  dimension  n  =  of 

the  square  imaging  grid. 

2.  Althou^  the  condition  nxiiiber  and,  hence,  the  error  bound  in  (46) 

are  exponentially  increasing  with  increasing  edge  dimension,  the  actual 
error  as  computed  by  (48)  does  not  show  this  increase  with  edge 
dimension. 

3.  The  amount  of  overdetermination  required  to  make  the  norm  error  in  the 
image  as  low  as  the  norm  error  in  the  scattered  field  is  extensive,  as 
much  as  300  percent. 

4.  For  weakly  overdetermined  systems,  increasing  the  detector  radius 
significantly  Increases  the  condition  ranter,  but  for  strongly  over¬ 
determined  systems,  it  makes  very  little  difference. 

5.  Increasing  the  percentage  overdetermination  improves  (decreases)  both 
the  ipper  bound  and  the  actual  norm  error  a  in  the  image. 


6.  For  a  fixed  radlijs  and  niniber  of  pixels  and  for  fixed  noise  in  the 
scattered  field,  the  condition  nxnher  and  noise  in  the  cooputed  image 
of  Y  both  decrease  as  the  contrast  in  the  test  object  is  increased. 

For  exanple,  for  detector  radius  R  =  10  and  300  percent  overdetermined 
data  and  for  a  7  by  7  image,  the  condition  nuiber  and  noise  are  860 
and  6  percent,  respectively,  for  a  test  object  consisting  of  a  centered 
Gaiissian  radial  distribution  having  a  peak  ccnplex  value  y  of  0.1- 
i  0.01,  but  the  condition  mnber  and  noise  are  200  and  4  percent, 
respectively,  for  the  same  test  object  with  Gaussian  peak  of  1.0-i  0.1  . 

SUMMARY 

Presently  clinical  ultrasomd  imaging  is  done  almost  exclusively  with  the 
B-scan  mode  (a  2-D  echo  strength  to  bri^tness  mode)  with  A-scan  (an  amplitude 
of  echo  graph  along  a  single  line)  or  with  Doppler  scanning  (for  measuring  blood 
flow  velocity) .  Ultrasomd  transmission  corputed  tomography  has  not  found 
clinical  application  yet  because  of  the  low  resolution  (no  better  than  5  to  10  mm) 
of  the  time  of  flight  and  attenuation  images.  An  improved  transmission  technique 
called  diffraction  tomograi^y  is  based  on  perturbation  solutions  to  the  Helmholtz 
wave  equation  and  in  principle  provide  improved  resolution  for  very  low  attenuating 
objects,  but  has  as  yet  not  been  implemented  due  to  significantly  large  attenuation 
in  tissues.  Diffraction  tomography  is  based  on  the  use  of  the  so  called  Bom  or 
i^ov  approximation  and  therefore  does  not  provide  a  complete  description  of  wave 
phenomena. 

We  have  developed  a  method  for  solving  an  exact  wave  integral  wave  equation 
for  an  inverse  scattering  solution  for  both  the  compressibility  and  density 
distribirtions  in  a  body.  Our  method  depends  on  the  \jse  of  multiple  incident  fields 
and  the  representation  of  the  total  fields  and  material  properties  in  a  basis 
function  e}q>ansion.  These  techniques  are  modifications  of  the  mcment  method  or 
Galerkin  method.  The  solutions  obtained  by  this  method,  vden  redundant  data  are 


used,  are  lailque  and  are  stable  even  with  noisy  data.  VJe  illustrate  the 
robustness  of  our  method  by  presenting  a  solution  image  made  with  10  percent 
noisy  data.  We  also  provide  tabulated  data  from  e:q)eriments  vdiich  show  the 
effect  of  picture  size  (niuiber  of  pixels)  and  degree  of  overdetermination 
on  the  image  quality.  The  present  algorithns  used  for  obtaining  a  solution 
require  an  enormous  amount  of  confutation  and  therefore  are  too  slow  and 
expensive  for  clinical  applications  [24] .  Inprovements  in  algorithms  speed  [24]  and 
anticipated  advances  in  computer  technology  could  open  applications  for 
laboratory  scanning  and  eventually  perhaps  even  for  clinical  scanning  for 
Galerkin  or  moment  method  solutions.  The  demonstration  of  the  accuracy  and 
robustness  of  these  new  methods  for  a  small  noiiber  of  pixels  should  act  as 
an  incentive  for  further  research  into  method  of  moment  or  other  more  exact 
approaches  to  inverse  scattering. 
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TABLE  1.  Upper  bound  on  the  norm  error  and  actual  norm  error  a  for 
various  size  images  as  function  of  percent  overdetermination  of 
the  problem.  Constant  grid  spacing  h  ■=  1/4  X  ,  problem  size  is 
NxN  ,  detector  radixis  is  at  10  wavelengths .  Five  percent  Gavissian 
noise  was  added  to  the  scattered  field.  A  two-dimensional  Gaussian 
test  object  similar  to  that  in  Figure  1  was  used. 
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FIGURE  1.  An  11  by  11  image  of  a  Gaussian  test  object.  The  pixels,  i.e.,  grid 
points,  are  separated  by  1/4  wavelength.  Thus,  the  image  is  11/4  by 
11/4  wavelengths.  Each  pair  of  upper  and  Icwer  nvirbers  corresponds 
to  a  scaled  valve  of  at  each  grid  point.  The  top  nijiber  in  each 
pair  is  10,000  times  the  value  of  the  real  part  of  and  the  bottom 
nvnber  is  10,000  times  the  imagiiiary  part  of  .  Using  this  test 
object,  simulated  scattering  data  was  obtained  using  the  method  of 
moments  to  solve  the  direct  solution.  Ten  percent  Gaussian  random  noise 
was  then  added  to  the  simulated  scattering  data,  i.e. ,  the  ratio  of  the 
one  norm  of  the  noise  to  the  one  norm  of  the  scattering  data  is  ten  per¬ 
cent.  The  simulated  scattering  data  with  noise  thus  generated  was  used 
as  input  for  the  method  of  moments  inverse  scattering  solution  method. 
The  results  of  the  calculations  to  retrive  the  original  object  are  shown 
in  the  next  figure. 
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FIGURE  2.  Reconstructed  image  of  object  shown  in  Figure  1  using  the  method 
of  moments  without  spatial  filtering.  The  image  was  made  frcm 
simulated  scattered  data  with  ten  percent  added  noise  using  the 
test  object  of  Figure  1  and  the  iterative  reconstruction  algorithm 
described  in  the  text.  The  above  results  after  20  iterations  shows 
about  17  percent  noise,  i.e. ,  a  «  0,17  from  equation  (A8). 
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FIGUEE  3.  Reconstxuction  of  image  shewn  in  Figure  1  using  the  method  of  moments 
and  spatial  filtering.  This  image  is  the  spatial  filtered  version  of 
Figure  2.  Since  Figure  2  is  sampled  at  k  sanples  per  wavelejigth  of 
the  incident  field  the  image  contains  noise  frequencies  tp  to  twice 
the  reciprocal  of  the  incident  fields  wave  length.  Let  Xq  be  the 
incident  field's  wavelength.  Then  spatial  freqxjencies  in  the  range 
from  1/Xq  to  2/Aq  may  be  ranoved  from  Figure  2  without  effecting 
the  practical  resolution  limit  of  the  image  of  .  Figure  3  is  the 
image  in  Figure  2  filtered  in  this  manner.  The  noise  is  reduced 
considerably  since  a  as  defined  by  equation  (48)  is  reduced  from 
0.17  in  Figure  2  to  0.05  here. 
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ABSTRACT 


This  paper  describes  a  Galerkin  scheme  v*iich  uses  a  sine  basis  (i.e., 
translations  of  sin(itx)/('nx)  —[15])  to  carry  out  an  ultrasonic  tomography 
inversion  based  on  the  eqiiation 

(1)  +  k^f  =  -  V.[Yp7f]  . 

In  (1) ,  f  denotes  the  sound  pressure, 

k  »  2Tr  K  frequency  x  =  (•c-*^)/k:q  ,  =  (p-Pq)/pq  vhere  ic 

(rsp.  Kq)  denotes  the  compressibility  of  the  body  (rsp.  the  fluid  surrounding 
the  body)  and  p  (rsp,  Pq)  denotes  the  corresponding  density.  The 
reconstruction  of  the  functions  y  and  y  in  1  (n  =  2  or  3)  is  based  on 

IC  p 

measuring  the  sound  pressure  f  on  a  finite  number  of  n-1  dimensional 
hyperplanes  for  each  fixed  source,  and  then  repeating  the  experiment  using  a 
finite  number  of  source  positions.  The  process  is  nonlinear,  since  all 
of  f  ,  Yp  and  Y^  Tnust  be  reconstructed  in  the  interior  of  the  body.  We 
thus  describe  an  iterative  procedure  for  carrying  out  the  reconstruction.  An 
example  is  described,  illustrating  a  two  dimensional  reconstruction  in  the 
presence  of  10%  Gaussian  noise. 
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